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Abstract 



Gaussian processes are rich distributions over functions, which provide a Bayesian nonparametric 
approach to smoothing and interpolation. We introduce simple closed form kernels that can be used 
with Gaussian processes to discover patterns and enable extrapolation. These kernels are derived by 
modelling a spectral density - the Fourier transform of a kernel - with a Gaussian mixture. The pro- 
posed kernels support a broad class of stationary covariances, but Gaussian process inference remains 
simple and analytic. We demonstrate the proposed kernels by discovering patterns and performing long 
range extrapolation on synthetic examples, as well as atmospheric CO2 trends and airline passenger 
data. We also show that we can reconstruct standard covariances within our framework. 



1 Introduction 

Machine learning is fundamentally about pattern discovery. The first machine learning models, such as 
the perceptron (Rosenblatt, 1962), were based on a simple model of a neuron (McCulloch and Pitts, 1943). 
Papers such as Rumelhart ct al. (1986) inspired hope that it would be possible to develop intelligent 
agents with models like neural networks, which could automatically discover hidden representations in 
data. Indeed, machine learning aims not only to equip humans with tools to analyze data, but to fully 
automate the learning and decision making process. 

Research on Gaussian processes (GPs) within the machine learning community developed out of neural 
networks research, triggered by Ncal (1996), who observed that Bayesian neural networks became Gaussian 
processes as the number of hidden units approached infinity. Neal (1996) conjectured that "there may be 
simpler ways to do inference in this case" . 

These simple inference techniques became the cornerstone of subsequent Gaussian process models for ma- 
chine learning (Rasmusscn and Williams, 2006). These models construct a prior directly over functions, 
rather than parameters. Assuming Gaussian noise, one can analytically infer a posterior distribution 
over these functions, given data. Gaussian process models have become popular for non-linear regres- 
sion and classification (Rasmusscn and Williams, 2006), and often have impressive empirical performance 
(Rasmusscn, 1996). 

The properties of likely functions under a GP, e.g., smoothness, periodicity, etc., are controlled by a 
positive definite covariance kernel^, an operator which determines the similarity between pairs of points in 
the domain of the random function. The choice of kernel profoundly affects the performance of a Gaussian 

*http : //mlg . eng. cam. ac .uk/andrew 
^http : / /people . seas .harvard. edu/ ~rpa/ 

^The terms covariance kernel, covariance function, kernel function, and kernel are used interchangeably. 
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process on a given task - as much as the choice of architecture, activation functions, learning rate, etc., 
can affect the performance of a neural network. 

Gaussian processes are sometimes used as expressive statistical tools, where the pattern discovery is per- 
formed by a human, and then hard coded into parametric kernels. Often, however, the squared exponential 
(Gaussian) kernel is used by default. In either case, Gaussian processes are used as smoothing interpola- 
tors with a fixed (albeit infinite) set of basis functions. Such simple smoothing devices are not a realistic 
replacement for neural networks, which were envisaged as intelligent agents that could discover hidden 
features in data^ via adaptive basis functions (MacKay, 1998). 

However, Bayesian nonparametrics can help build automated intelligent systems that reason and make 
decisions. It has been suggested that the human ability for inductive reasoning - concept generalization 
with remarkably few examples - could derive from a prior combined with Bayesian inference (Yuillc and 
Kersten, 2006; Tencnbaum et al., 2011; Stcyvcrs ct al., 2006). Bayesian nonparametric models, and Gaus- 
sian processes in particular, are an expressive way to encode prior knowledge, and also reflect the belief 
that the real world is infinitely complex (Neal, 1996). 

With more expressive kernels, one could use Gaussian processes to learn hidden representations in data. 
Expressive kernels have been developed by combining Gaussian processes in a type of Bayesian neural 
network structure (Salakhutdinov and Hinton, 2008; Wilson et al., 2012; Damianou and Lawrence, 2012). 
However, these approaches, while promising, typically 1) are designed to model specific types of structure 
(e.g., input-dependent correlations between different tasks); 2) make use of component GPs with simple 
interpolating kernels; 3) indirectly induce complicated kernels that do not have a closed form and are 
difficult to interpret; and 4) require sophisticated approximate inference techniques that are much more 
demanding than that required by simple analytic kernels. 

Sophisticated kernels are most often achieved by composing together a few standard kernel functions (Ar- 
chambeau and Bach, 2011; Durrandc ct al., 2011; Gonen and Alpaydm, 2011; Rasmussen and Williams, 
2006). Tight restrictions are typically enforced on these compositions and they are hand-crafted for spe- 
cialized applications. Without such restrictions, complicated compositions of kernels can lead to overfitting 
and unmanageable hyperparameter inference. Moreover, while some compositions (e.g., addition) have an 
easily interpretable effect on the resulting sample paths, many other operations change the distribution 
over functions in ways that are difficult to identify. It is difficult, therefore, to construct an effective induc- 
tive bias for kernel composition that leads to automatic discovery of the appropriate statistical structure, 
without human intervention. 

This difficulty is exacerbated by the fact that it is challenging to say anything about the covariance function 
of a stochastic process from a single draw if no assumptions are made. If we allow the covariance between 
any two points in the input space to arise from any positive definite function, then we gain essentially 
no information from a single realization. Most commonly one assumes a restriction to stationary kernels, 
meaning that covariances are invariant to translations in the input space. 

In this paper, we explore flexible classes of kernels that go beyond composition of simple analytic forms, 
while maintaining the useful inductive bias of stationarity. We propose new kernels which can be used 
to automatically discover patterns and extrapolate far beyond the available data. This class of kernels 
contains many stationary kernels, but has a simple closed form that leads to straightforward analytic 
inference. The simplicity of these kernels is one of their strongest qualities. In many cases, these kernels can 
be used as a drop in replacement for the popular squared exponential kernel, with benefits in performance 
and expressiveness. By learning features in data, we not only improve predictive performance, but we 
can more deeply understand the structure of the problem at hand - greenhouse gases, air travel, heart 
physiology, brain activity, etc. 

After a brief review of Gaussian processes in Section 2, we derive the new kernels in Section 3 by modelling 
a spectral density with a mixture of Gaussians. We focus our experiments in Section 4 on elucidating 

^In this paper, wc refer to representations, features and patterns interchangeably. In other contexts, the term features 
sometimes specifically means low dimensional representations of data, like neurons in a neural network. 
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the fundamental differences between the proposed kernels and the popular alternatives in Rasmussen and 
Williams (2006). In particular, we show how the proposed kernels can automatically discover patterns and 
extrapolate on the CO2 dataset in Rasmussen and Williams (2006), on a synthetic dataset with strong 
negative covariances, on a difficult synthetic sine pattern, and on airline passenger data. We also use our 
framework to reconstruct several popular standard kernels. 

2 Gaussian Processes 

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian 
distribution. Using a Gaussian process, we can define a distribution over functions f{x), 

f(x)^gV[m{x),k{x,x')), (1) 

where x G is an arbitrary input variable, and the mean function m{x) and covariance kernel k{x,x') 
are defined as 

mix) = E[/(x)] , (2) 

k{x,x')^cov{f{x)J{x')). (3) 

Any collection of function values has a joint Gaussian distribution 

[f{x,)J{x2),...JixN)V ^Mifi,K), (4) 

where the N x N covariance matrix K has entries K^j = k{xi, Xj), and the mean fi has entries fii — m{xi). 
The properties of the functions - smoothness, periodicity, etc. - are determined by the kernel function. 

The popular squared exponential (SE) kernel has the form 

ksE{x,x') =exp(-0.5||a;-a;'||V^2)^ (5) 

Functions drawn from a Gaussian process with this kernel function are infinitely differentiable, and can 
display long range trends. Gaussian processes with a squared exponential kernel are simply smoothing 
devices: the only covariance structure that can be learned from data is the length-scale £, which determines 
how quickly a Gaussian process function varies with x. 

Assuming Gaussian noise, one can analytically infer a posterior predictive distribution over Gaussian 
process functions, and analytically derive a marginal likelihood of the observed function values y given 
only hyperparameters 6, and the input locations {xn}n=ij p{y\(^i {xn}n=i)- This marginal likelihood can 
be optimised to estimate hyperparameters such as £, or used to integrate out the hyperparameters via 
Markov chain Monte Carlo (Murray and Adams, 2010). Detailed Gaussian process references include 
Rasmussen and WiUiams (2006), Stein (1999), and Cressie (1993). 

3 Kernels for Pattern Discovery 

In this section we introduce a class of kernels that can discover patterns, extrapolate, and model negative 
covariances. This class contains a large set of stationary kernels. Roughly, a kernel measures the similarity 
between data points. As in Equation (3), the covariance kernel of a GP determines how the associated 
random functions will tend to vary with inputs (predictors) x € M^. A stationary kernel is a function of 
T = x ^ x' , i.e., it is invariant to translation of the inputs. 

Any stationary kernel (aka covariance function) can be expressed as an integral using Bochner's theorem 
(Bochncr, 1959; Stein, 1999; Gikhman and Skorokhod, 2004): 
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Theorem 3.1 (Bochner) A complex-valued function k on is the covariance function of a weakly sta- 
tionary mean square continuous complex-valued random process on M.^ if and only if it can be represented 
as 

k{T) = / e2--"-V(ds) , (6) 

where ^ is a positive finite measure. 

If tp has a density S'(s), then S is called the spectral density or power spectrum of fc, and k and S are 
Fourier duals (Chatficld, 1989): 



k{T) = J 5(s)e^"" ^ds , (7) 
S{s) = / k{T)e-^''''^^dT . (8) 



In other words, a spectral density entirely determines the properties of a stationary kernel. Substituting the 
squared exponential kernel of (5) into (8), we find its spectral density is S'se(s) = (27r£^)^/^ exp(— 27r^£^s^). 
Therefore SE kernels, and mixtures of SE kernels, are a very small corner of the set of possible stationary 
kernels, as they correspond only to Gaussian spectral densities centered on the origin. 

However, by using a mixture of Gaussians that have non-zero means, one can achieve a much wider range 
of spectral densities. Indeed, mixtures of Gaussians are dense in the set of all distribution functions 
(Kostantinos, 2000). Therefore, the dual of this set is also dense in stationary covariances. That is, we can 
approximate any stationary covariance kernel to arbitrary precision, given enough mixture components in 
the spectral representation. This observation motivates our approach, which is to model GP covariance 
functions via spectral densities that are scale-location mixtures of Gaussians. 

We first consider a simple case, where 

0(s;A^,cr^)= J_ exp{--'^(s-pt)^}, and (9) 

S{s)^ms) + cf>{-s)]/2, (10) 

noting that spectral densities are real and symmetric about s — (since A:(r) is real and symmetric about 
T = 0) (e.g. Hormander (1990)). 

Substituting S{s) into equation (7), we find 

fc(r) = exp{-27r2T^cr^} cos(27rr/i) . (11) 

If (/)(s) is instead a mixture of Q Gaussians on M^, where the g'^ component has mean vector fiq — {nq^\ . . . , Mg^'') 
and covariance matrix Mg = diag{vq^\ . . . ,Vq^^), and Tp is the p*'' component of the P dimensional vec- 
tor T ^ X ~ x' , then 

Q P 

fc(r) = ^w;,n exp{-27rV>(^')}cos(27rrpAi(f)). (12) 

The integral in (7) is tractable even when the spectral density is an arbitrary Gaussian mixture, allowing 
us to derive'' the exact closed form expressions in Eqs. (11) and (12), and to perform analytic inference 
with Gaussian processes. Moreover, this class of kernels is expressive - containing many stationary kernels 
- but nevertheless has a simple form. 

These kernels are easy to interpret, and provide drop-in replacements for kernels in Rasmussen and 
Williams (2006). The weights Wq specify the relative contribution of each mixture component. The 



^Detailed derivations of Eqs. (11) and (12) are in the Appendix. 
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inverse means 1/nq are the component periods, and the inverse standard deviations are length- 

scales, determining how quickly a component varies with the inputs x. The kernel in Eq. (12) can also be 
interpreted through its associated spectral density. In Section 4, we use the learned spectral density to 
interpret the number of discovered patterns in the data and how these patterns generalize. Henceforth, 
we refer to the kernel in Eq. (12) as a spectral mixture (SM) kernel. 



4 Experiments 

We show how the SM kernel in Eq. (12) can be used to discover patterns, extrapolate, and model negative 
covariances. We contrast the SM kernel with popular kernels in, e.g., Rasmussen and Williams (2006) and 
Abrahamscn (1997), which typically only provide smooth interpolation. Although the SM kernel generally 
improves predictive likelihood over popular alternatives, we focus on clearly visualizing the learned kernels 
and spectral densities, examining patterns and predictions, and discovering structure in the data. Our 
objective is to elucidate the fundamental differences between the proposed SM kernel and the alternatives. 

In all experiments, Gaussian noise is assumed, so that marginalization over the unknown function can 
be performed in closed form. Kernel hyperparameters are trained using nonlinear conjugate gradients to 
optimize the marginal likelihood p{y\0,{xn}n=i) of the data y given hyperparameters 9, as described in 
Section 2, assuming a zero mean GP. Automatic relevance determination (MacKay et al., 1994; Tipping, 
2004) takes place during training, removing extraneous components from the proposed model, through 
the complexity penalty in the marginal likelihood (Rasmussen and Williams, 2006). For a fully Bayesian 
treatment, the spectral density could alternatively be integrated out using Markov chain Monte Carlo 
(Murray and Adams, 2010), rather than choosing a point estimate. However, we wish to emphasize that 
the SM kernel can be successfully used in the same way as other popular kernels, without additional 
inference efforts. 

We compare with the popular squared exponential (SE), Matern (MA), rational quadratic (RQ), and 
periodic (PE) kernels. In each comparison, we attempt to give these alternative kernels fair treatment: 
we initialise hyperparameters at values that give high marginal likelihoods and which are well suited to 
the datasets, based on features we can already see in the data. Conversely, we randomly initialise the 
parameters for the SM kernel. Training runtimes are on the order of minutes for all tested kernels. 



4.1 Extrapolating Atmospheric CO2 

Keeling and Whorf (2004) recorded monthly average atmospheric CO2 concentrations at the Mauna Loa 
Observatory, Hawaii. The data are shown in Figure 1. The first 200 months are used for training (in blue), 
and the remaining 301 months (~ 25 years) are used for testing (in green). 

This dataset was used in Rasmussen and Williams (2006), and is frequently used in Gaussian process 
tutorials, to show how GPs are flexible statistical tools: a human can look at the data, recognize patterns, 
and then hard code those patterns into covariance kernels. Rasmussen and Williams (2006) identify, by 
looking at the blue and green curves in Figure la, a long term rising trend, seasonal variation with possible 
decay away from periodicity, medium term irregularities, and noise, and hard code a stationary covariance 
kernel to represent each of these features. 

However, in this view of GP modelling, all of the interesting pattern discovery is done by the human 
user, and the Gaussian process is used simply as a smoothing device, with the flexibility to incorporate 
human inferences in the prior. Our contention is that such pattern recognition can also be performed 
algorithmically. To discover these patterns without encoding them a priori into the GP, we use the 
spectral mixture kernel in Eq. (12), with Q = 10. The results are shown in Figure la. 
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Figure 1: Mauna Loa CO2 Concentrations, a) Forecasting CO2. The training data are in blue, the testing data in 
green. Mean forecasts made using the SM kernel are in black, with 2 stdev about the mean (95% of the predictive 
mass) in gray shade. Predictions using the Matern (MA), squared exponential (SE), rational quadratic (RQ), and 
periodic kernels (PE) are in cyan, dashed red, magenta, and orange, respectively, b) The log spectral densities of 
the learned SM and squared exponential kernels are in black and red, respectively. 



In the training region, predictions using each kernel are essentially equivalent, and entirely overlap with 
the training data. However, unlike the other kernels, the SM kernel (in black) is able to discover patterns 
in the training data and accurately extrapolate over a long range. The 95% high predictive density (HPD) 
region contains the true CO2 measurements for the duration of the measurements. 

We can see the structure discovered by the SM kernel in the learned log spectral density in Figure lb. Of 
the Q — 10 components, only seven were used. This is an example of automatic relevance determination 
(ARD) occurring through the mixture weights. There is a peak at a frequency near 0.07, which corresponds 
to a period of about 14 months, roughly a year, which is in line with seasonal variations. However, there 
are departures away from periodicity, and accordingly, there is a large variance about this peak, reflecting 
the uncertainty in this learned feature. We see other, sharper peaks corresponding to periods of 6 months, 
4 months, 3 months and 1 month. 

In red, we show the spectral density for the learned SE kernel, which covers some of the peaks in the SM 
spectral density, but assigns significant mass to many frequencies which are not important for explaining 
the data. All SE kernels have a Gaussian spectral density centred on zero. Since a mixture of Gaussians 
centred on zero is a poor approximation to many more general density functions, combinations of SE 
kernels have limited expressive power. Indeed the spectral density learned by the SM kernel in Figure lb 
is highly non-Gaussian. The test predictive performance using the SM, SE, MA, RQ, and PE kernels is 
given in Table 1, under the heading CO2. 



4.2 Recovering Popular Kernels 

The SM class of kernels contains many stationary kernels, since mixtures of Gaussians can be used to 
construct a wide range of spectral densities. Even with a small number of mixture components, e.g., Q < 10, 
the SM kernel can closely recover the most popular stationary kernels catalogued in Rasmussen and 
WiUiams (2006). 

As an example, we start by sampling 100 points from a one-dimensional GP with a Matern kernel with 
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Figure 2: Recovering popular correlation functions (normalised kernels), with t = x — x' . The true correlation 
function underlying the data is in green, and SM, SE, and empirical correlation functions are in dashed black, red, 
and magenta, respectively. Data are generated from a) a Matern kernel, and b) a sum of rational quadratic and 
periodic kernels. 



degrees of freedom = 3/2: 




where ^ = 5 and a = 4. Sample functions from a Gaussian process with this Matern kernel are far less 
smooth (only once-differentiable) than Gaussian process functions with a squared exponential kernel. 

We attempt to reconstruct the kernel underlying the data by training an SM kernel with Q — 10. After 
ARD in training, only (5 = 4 components are used. The log marginal likelihood of the data - having 
integrated away the Gaussian process - using the trained SM kernel is —133, compared to —138 for the 
Matern kernel that generated the data. Training the squared exponential kernel in (5) gives a log marginal 
likehhood of -140. 

Figure 2a shows the learned SM correlation function"*, compared to the generating Matern correlation 
function, the empirical autocorrelation function, and learned squared exponential correlation function. 
Although often used in geostatistics to guide choices of Gaussian process kernel functions (and parameters) 
(Cressie, 1993), the empirical autocorrelation function tends to be unreliable, particularly with a small 
amount of data (e.g., N < 1000), and at high lags (for r > 10). In Figure 2a, the empirical autocorrelation 
function is erratic and does not resemble the Matern kernel for r > 10. Moreover, the squared exponential 
kernel cannot capture the heavy tails of the Matern kernel, no matter what length-scale it has. The learned 
SM kernel more closely matches the Matern kernel. 

Next, we reconstruct a mixture of the rational quadratic (RQ) and periodic kernels (PE) in Rasmusscn 
and Williams (2006): 

ARQ(r) = (l + -^)-", (14) 
fcpE(T) = exp(-2sin2(7rra;)/£|,£;) . (15) 

The rational quadratic kernel in (14) is derived as a scale mixture of squared exponential kernels with 
different length-scales. The standard periodic kernel in (15) is derived by mapping the two dimensional 

*A correlation function c(x,x') is a normalised covariance kernel k{x,x'), such that c{x,x') = k(x,x'^ j ^k(x^ x)k{x', x') 
and c{x, x) = 1. 
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Figure 3: Negative Covariances. a) Observations of a discrete time autoregressive (AR) series with negative 
covariances. b) The SM learned kernel is shown in black, while the true kernel of the AR series is in green, with 
T = X — x' . c) The spectral density of the learned SM kernel is in black. 

variable u{x) = (cos(x), sin(x)) through the squared exponential kernel in (5). Derivations for both the RQ 
and PE kernels in (14) and (15) are in Rasmussen and Williams (2006). Rational quadratic and Matern 
kernels are also discussed in Abrahamsen (1997). We sample 100 points from a Gaussian process with 
kernel lOfcRQ + 4A:pe, with a = 2, w = 1/20, (.rq = 40, ^PE = 1- 

We reconstruct the kernel function of the sampled process using an SM kernel with Q = 4, with the 
results shown in Figure 2b. The heavy tails of the RQ kernel are modelled by two components with large 
periods (essentially aperiodic): one of these components has a short length-scale, while the other has a 
large length-scale. The third component has a relatively short length-scale, and a period of 20. There is 
not enough information in the 100 sample points to justify using more than (5 = 3 components, and so - 
through automatic relevance determination in training - the fourth component in the SM kernel moves 
towards zero. The empirical autocorrelation function somewhat captures the periodicity in the data, but 
significantly underestimates the correlations. The squared exponential kernel learns a long length-scale: 
since the SE kernel is significantly misspccificd with the true generating kernel, the data are explained as 
noise. 



4.3 Negative Covariances 

All of the stationary covariance functions in the standard machine learning Gaussian process reference Ras- 
mussen and Williams (2006) are everywhere positive, including the periodic kernel, A;(r) = exp(— 2 sin^(7r r w)/£^). 
While positive covariances are often suitable for interpolation, capturing negative covariances can be es- 
sential for extrapolating patterns: for example, linear trends have long-range negative covariances. We 
test the ability of the SM kernel to learn a negative covariance function, by sampling 400 points from a 
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simple AR(1) discrete time Gaussian process: 



y{x + l) 



— e 



0.01 



y{x) + (T€{x) , 



(16) 
(17) 



e(a;) ^ AA(0, 1) , 



which has kernel 



k{x, x) 




'l)k-x'l/(i_g-.02) 



(18) 



The process in Eq. (16) is shown in Figure 3a. This process follows an oscillatory pattern, systematically 
switching states every x = 1 unit, but is not periodic and has long range covariances: if we were to only 
view every second data point, the resulting process would vary rather slowly and smoothly. 

We see in Figure 3b that the learned SM covariance function accurately reconstructs the true covariance 
function. The spectral density in Figure 3c) shows a sharp peak at a frequency of 0.5, or a period of 2. 
This feature represents the tendency for the process to oscillate from positive to negative every x — \ unit. 
In this case (3 = 4, but after automatic relevance determination in training, only 3 components were used. 
Using the SM kernel, we forecast 20 units ahead and compare to other kernels in Table 1 (NEG CDV). 



4.4 Discovering the Sine Pattern 

The sine function is defined as sinc(x) = sin(7ra;)/(7ra;). We created a pattern by combining three sine 
functions: 



This is a complex oscillatory pattern. Given only the points shown in Figure 4a, we wish to complete the 
pattern for x G [—4.5, 4.5]. Unlike the CO2 example in Section 4.1, it is perhaps even difficult for a human 
to extrapolate the missing pattern from the training data. It is an interesting exercise to focus on this 
figure, identify features, and fill in the missing part. 

Notice that there is complete symmetry about the origin a: = 0, peaks at a; = — 10 and x = 10, and 
destructive interference on each side of the peaks facing the origin. We therefore might expect a peak 
at x = and a symmetric pattern around x — 0. 

As shown in Figure 4b, the SM kernel with Q = 10 reconstructs the pattern in the region x E [—4.5,4.5] 
almost perfectly from the 700 training points in blue. Moreover, using the SM kernel, 95% of the posterior 
predictive mass entirely contains the true pattern in x G [—4.5,4.5]. GPs using Matern, squared expo- 
nential, rational quadratic, and periodic kernels are able to predict reasonably within x = 0.5 units of the 
training data, but entirely miss the pattern in x e [—4.5, 4.5]. 

Figure 4c shows the learned SM correlation function (normalised kernel). For t e [0, 10] there is a local 
pattern, roughly representing the behaviour of a single sine function. For r > 10 there is a repetitive 
pattern representing a new sine function every 10 units - an extrapolation a human might make. With a 
sine functions centred at x ~ —10, 0, 10, we might expect more sine patterns every 10 units. The learned 
Matern correlation function is shown in red in Figure 4c - unable to discover complex patterns in the data, 
it simply assigns high correlation to nearby points. 

Figure 4d shows the (highly non-Gaussian) spectral density of the SM kernel, with peaks at 0.003, 0.1, 
0.2, 0.3, 0.415, 0.424, and 0.492. In this case, only 7 of the Q = 10 components are used. The peak 
at 0.1 represents a period of 10: every 10 units, a sine function is repeated. The variance of this peak is 
small, meaning the method will extrapolate this structure over long distances. By contrast, the squared 
exponential spectral density simply has a broad peak, centred at the origin. The predictive performance 
for recovering the missing 300 test points (in green) is given in Table 1 (SING). 



y{x) = sinc(x + 10) + sinc(a;) + sinc(a; — 10) . 



(19) 
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Figure 4: Discovering a sine pattern, a) The training data are shown in blue. The goal is to fill in the missing 
region x £ [—4.5,4.5]. b) The training data are in blue, the testing data in green. The mean of the predictive 
distribution using the SM kernel is shown in dashed black. The mean of the predictive distributions using the 
squared exponential, Matern, rational quadratic, and periodic kernels, are respectively in red, cyan, magenta, and 
orange, c) The learned SM correlation function (normalised kernel) is shown in black, and the learned Matern 
correlation function is in red, with t = x — x' . d) The log spectral densities of the SM and squared exponential 
kernels are in black and red, respectively. 
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Figure 5: Predicting airline passenger numbers, a) The training and testing data are in blue and green, respectively. 
The mean of the predictive distribution using the SM kernel is shown in black, with 2 standard deviations about the 
mean (95% of the predictive probability mass) shown in gray shade. The mean of the predictive distribution using 
the SM kernel is shown in black. The mean of the predictive distributions using the squared exponential, Matern, 
rational quadratic, and periodic kernels, are in red, cyan, magenta, and orange, respectively. In the training region, 
the SM and Matern kernels are not shown, since their predictions essentially overlap with the training data, b) 
The log spectral densities of the SM and squared exponential kernels are in black and red, respectively. 

4.5 Airline Passenger Data 

Figure 5a shows airline passenger numbers, recorded monthly, from 1949 to 1961 (Hyndman, 2005). Based 
on only the first 96 monthly measurements, in blue, we wish to forecast airline passenger numbers for the 
next 48 months (4 years); the corresponding 48 test measurements are in green. Companies wish to make 
such long range forecasts to decide upon expensive long term investments, such as large passenger jets, 
which can cost hundreds of millions of dollars. 

There are several features apparent in these data: short seasonal variations, a long term rising trend, and 
an absence of white noise artifacts. Many popular kernels are forced to make one of two choices: 1) Model 
the short term variations and ignore the long term trend, at the expense of extrapolation. 2) Model the 
long term trend and treat the shorter variations as noise, at the expense of interpolation. 

As seen in Figure 5a, the Matern kernel is more inclined to model the short term trends than the smoother 
SE or RQ kernels, resulting in sensible interpolation (predicting almost identical values to the training 
data in the training region), but poor extrapolation - moving quickly to the prior mean, having learned 
no structure to generalise beyond the data. The SE kernel interpolates somewhat sensibly, but appears to 
underestimate the magnitudes of peaks and troughs, and treats repetitive patterns in the data as noise. 
Extrapolation using the SE kernel is poor. The RQ kernel, which is a scale mixture of SE kernels, is more 
able to manage different length-scales in the data, and generalizes the long term trend better than the SE 
kernel, but interpolates poorly. 

By contrast, the SM kernel interpolates nicely (overlapping with the data in the training region), and is 
able to successfully extrapolate complex patterns far beyond the data, capturing the true airline passenger 
numbers for years after the data ends, within a small band containing 95% of the predictive probability 
mass. 

Of the Q = 10 initial components specified for the SM kernel, 7 were used after ARD in training. The 
learned spectral density in Figure 5b shows a large sharp low frequency peak (at about 0.00148). This 
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Table 1: We compare the test performance of the proposed spectral mixture (SM) kernel with squared 
exponential (SE), Matern (MA), rational quadratic (RQ), and periodic (PE) kernels. The SM kernel 
consistently has the lowest mean squared error (MSE) and highest log likelihood (£). 





SM 


SE 


MA 


RQ 


PE 


C02 


MSE 


9.5 


1200 


1200 


980 


1200 




170 


-320 


-240 


-100 


-1800 


NEG GOV 


MSE 


62 


210 


210 


210 


210 


C 


-25 


-70 


-70 


-70 


-70 


SING 


MSE 


0.000045 


0.16 


0.10 


0.11 


0.05 


C 


3900 


2000 


1600 


2000 


600 


AIRLINE 


MSE 


460 


43000 


37000 


4200 


46000 


C 


-190 


-260 


-240 


-280 


-370 



peak corresponds to the rising trend, which generalises well beyond the data (small variance peak), is 
smooth (low frequency), and is important for describing the data (large relative weighting). The next 
largest peak corresponds to the yearly trend of 12 months, which again generalises, but not to the extent 
of the smooth rising trend, since the variance of this peak is larger than for the peak near the origin. The 
higher frequency peak &t x = 0.34 (period of 3 months) corresponds to the beginnings of new seasons, 
which can explain the effect of seasonal holidays on air traffic. We could interpret the peak at x = 0.25 
as relating to seasonal holidays as well, since there is a major holiday every 4 months. A more detailed 
study of these features and their properties - frequencies, variances, weightings, etc. - could isolate less 
obvious non-seasonal effects important for modelling airline passenger numbers. Table 1 (AIRLINE) shows 
comparative predictive performance for forecasting 48 months of airline passenger numbers. 



5 Discussion 

We have derived expressive closed form kernels and we have shown that these kernels, when used with 
Gaussian processes, can discover patterns in data and extrapolate over long ranges. The simplicity of 
these kernels is one of their strongest properties: they can be used as drop-in replacements for popular 
kernels such as the squared exponential kernel, with major benefits in expressiveness and performance, 
while retaining simple training and inference procedures. 

Gaussian processes have proven themselves as powerful smoothing interpolators. We believe that pattern 
discovery and extrapolation is an exciting new direction for Bayesian nonparametric approaches, which can 
capture rich variations in data. Here we have shown how Bayesian nonparametric models can naturally 
be used to generalise a pattern from a small number of examples. 

We have only begun to explore what could be done with such pattern discovery methods. In future work, 
one could integrate away a spectral density, using recently developed efficient MCMC for GP hyperparame- 
ters (Murray and Adams, 2010). Moreover, spectral representations of the proposed models suggest recent 
Toeplitz methods (Saatchi, 2011) could be applied for significant speedup in inference and predictions. 
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6 Appendix 



6.1 Derivations for Spectral Mixture Kernels 



A stationary kernel k{x,x') is the inverse Fourier transform of its spectral density S{s), 




(20) 



where t = x — x' . First suppose 



S{s) = 



1 




(21) 
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where s, fi, a and t = x — x' are scalars. Substituting (21) into (20), 

k(x,x')= exp(2TTis(x - x')) , exp(-—r(s - ^)'^)ds (22) 
J V27rCT2 2cr^ 



let T — X — x' 



If 1 

exp[27rir - ^ (s^ - 2/is + fj.'^)]ds (23) 



V27rCT2 y 2(7 



1 /" r 1 2 , /n ■ , M X A*' 



2 



exp[-— + (27riT+ -^)s" ^]ds (24) 



J 2cr 

let a = -—7 , 27rir H — - , c = - -— r 

2(T-^ CT^ 2(7^ 

= 72b/'"p^""^'^"^^'^'"p^£+^^'' ^''^ 

^exp[(2vr^r+i^)^^^-i^] (26) 

= exp[(~47rV2 + 47r*ri^ + g)^ _ i^] (27) 
= exp[-27r2(a; - x')^cr^] [cos(27r(a; - x')^) + i sin(27r(x - x')fi))] . (28) 
Noting that the spectral density S{s) must be symmetric about s — 0, we let 

<l){s;fi,(T^) = ^'_ exp{--j^(s-pt)^}, and (29) 

S{s) = [4>is) + cf>{-s)]/2. (30) 
Closely following the above derivation, substituting (30) into (20) gives 

k{T) = exp{-27r W^} cos{2ttth) . (31) 



If (/)(s) is instead a mixture of Q Gaussians on M^, where the g"^ component has mean vector /Xg — {nq^\ ■ ■ ■ , Mg^^) 
and covariance matrix Mg = diag(wq^'', . . . ,Wg^-'), and Tp is the p*'' component of the P dimensional vec- 
tor T = X — x' , then the integral in (20) becomes a sum of a product of the one dimensional integrals we 
encountered to derive (31), from which it follows that 

Q P 

fc(r) = J2'^,l[eM-^^'ryP^}cos{2nrpfiil^^). (32) 

9=1 P=l 

6.2 Comment on Training Hyperparameters 

Generally, we have had success naively training kernel hyperparameters using conjugate gradients (we use 
Carl Rasmussen's 2010 version of minimize . m) to maximize the marginal likelihood p{y\0) of the data 
y given hypers 9, having analytically integrated away a zero mean Gaussian process. We have found 
subtracting an empirical mean from the data prior to training hyperparameters (with conjugate gradients) 
undesirable, sometimes leading to local optima with lower marginal likelihoods, particularly on small 
datasets with a rising trend. 



